1. Field of the Invention
The invention relates to a method for forward modeling the spatial distribution faults and fractures in a geologic formation.
2. Background of the Invention
The prediction of faulting and fracturing is very important in oil and gas exploration and production. Seismic data is often used to find faults that bound or delineate hydrocarbon reservoirs. However, due to the limited resolution of seismic data, the details of the faulting in the subsurface may not be determinable.
Knowledge of the distribution of the fractures in a geologic formation is of great importance first for optimizing the location and the spacing between the wells that are to be drilled through an oil formation. Furthermore, the geometry of the fracture network conditions influences the displacement of fluids on the reservoir scale as well as on the local scale, where it determines the elementary matrix blocks in which the oil is trapped. Knowledge of the distribution of the fractures is thus also very useful at a later stage for the reservoir engineer who wants to extrapolate the production curves and to calibrate the models simulating reservoirs. The development of naturally fractured reservoirs thus requires better knowledge of the geometry of the fracture networks and of their contribution to the orientation of the flows.
Seismic data are commonly used for acquiring information about the subsurface of the Earth. Changes in the elastic properties of subsurface rocks appear as seismic reflections. Such changes in the properties of the rocks typically occur at boundaries between geologic formations, at fractures and at faults. The vertical resolution of the seismic method is approximately one-quarter wavelength of the seismic wave and, in typical situations, is of the order of 10 meters. The horizontal resolution is determined by the size of the Fresnel zone for the seismic wave at the depth of interest and may be tens or even hundreds of meters. By using sophisticated processing techniques, such as prestack migration taking advantage of data redundancy, the positions of the seismic reflectors may be more accurately determined up to the spacing of the geophones. Nevertheless, it is only the major seismic reflectors that may be determined by this method. Additional information about finer scale faulting and fracturing would be very useful in predicting flow characteristics of a hydrocarbon reservoir and in development of a program for optimization of hydrocarbon recovery. Since all faults are not clearly delineated by the seismic, it would also be useful to verify the location of faults that are difficult to see seismically are interpreted in a consistent manner.
U.S. Pat. No. 5,953,680 issued to Divies et al describes a method for creating a two-dimensional (2-D) kinematic model of a geologic basin affected by faults. The basin is divided into a number of layers or banks whose geometric positions are known. The tectonic deformation of each modeled layer is determined separately by taking its thickness and length into account, with compaction being taken into account. The basic assumption is that the banks are competent units that undergo little deformation. The method does not include the material properties of the rocks as part of the input and hence is not particularly well suited for determining the effects of loading.
U.S. Pat. No. 5,838,634 issued to Jones et al obtains geologic models of the subsurface that are optimized to match as closely as feasible geologic constraints known or derived from observed geologic data. The models also conform to geophysically based constraints indicated by seismic survey data. It accounts for geophysical information by converting the geologic model to synthetic seismic traces, accounting for fluid saturation, and comparing these traces with observed seismic trace data. The process perturbs the rock properties in the geologic model until the geologic model is consistent with geologic and geophysical data and interpretations. However, the issue of how to obtain a reasonable fine-scale geologic model is not addressed.
Broadly spealdng, four different ways have been used to model the problem of geologic modeling of the subsurface on a wide range of scales. The first approach is to use statistics to capture the number and orientation of faults at one scale or in one deformational setting and to use simple statistical rules to extrapolate this information to other scales or deformational settings. An example of this is U.S. Pat. No. 5,659,135 to Cacas.
A second method that has been used is to use finite element modeling to solve the stress field from given input deformations. Once stress exceeds a given amount a fault or fracture is drawn in by hand and then the model simulation can continue. Alternatively, faulting patterns are put in by hand, and the formation is pressured up to estimate a stress distribution. The rock is modeled is a network of distinct elastic elements, connected by elastic connection to its outer boundaries. The main obstacles to the application of such methods for geologic modeling are the computer time and the human interaction that is involved. The computer time roughly increases as the square of the number of nodes in the model and the models must be continuously interacted with by the user to put in new faults as they are believed to have occurred.
In a third method, large scale rules of geometry or faulting seen in the subsurface under certain deformation conditions are quantified and applied to forward modeling software. These forward models usually consist of a well-defined set of large scale shapes that are expected to be produced. An example of this is U.S. Pat. No. 5,661,698 issued to Cacas, which starts out with a group of major faults detected by means of an exploration of the zone, and additional minor faults that have not been detected during the exploration. The fractal characteristics of the major faults are determined and the additional minor faults are constrained to have the same fractal characteristics. The fractal characteristics used include the fractal dimension of the fault network and a density function defining a distribution of lengths of the faults. Such a method does not account for differences in the rock properties of different geologic formations and differences in their mode of faulting.
A fourth method that has been used is the so-called xe2x80x9cdistinct element model.xe2x80x9d It uses small scale rules of stress and strain to move nodes in a model to predict faulting and fracturing. It is well suited for problems of geologic fracturing but suffers from the drawback of being computationally slow. In addition, the methods are not particularly user friendly in terms of user interface used for specification of the model and of the material properties.
There is a need for an invention that is able to simulate faulting and fracturing on a variety of scales in a subsurface geological model. Such an invention should preferably take into account the differences in material strength of different types of rocks. Such an invention should also be computationally fast. In addition, it is preferable that the invention should be user friendly in that specification of the rock properties and loading be easily input and that the invention be able to provide graphical displaces of the deformation process. The present invention satisfies this need.
The present invention is a computer-implemented method for modeling of faulting and fracturing that uses xe2x80x9csmall scale rulesxe2x80x9d to produce large scale results. Organizationally, the software for the invention is made up of two parts. One part is a user interface for inputting deformations, pre-existing faults and fractures, and material rock properties. The second part of the software is the code that solves the motion of each point or node in the subsurface volume defined by the user interface. The solution code solves the forces on each node and their resulting movement. Faults and fractures occur as the nodes are widely separated and the forces between the nodes (that are based on the node locations) are changed. The user interface may be used to produce a quick look at the deformation results and to view the results of the full simulation.
Functionally, as with all modeling programs, the present invention has two primary components. The first is the model definition and the second is the response of the model to the applied deformation and stresses. In the present invention, there is a third important component: a graphical user interface that is very useful in both the model definition and in viewing the results of the applied deformation and stresses.
The model may be defined in one of three modes. In the aerial mode, the model is 2-dimensional with the material and a substrate on a horizontal plane. A nodal pattern that, in the preferred embodiment, is a regular triangular lattice is used for the model definition, although other nodal patterns may also be used. An alternate embodiment of the invention allows for randomly positioned nodes in a so-called xe2x80x9cdisordered lattice.xe2x80x9d The deformation can be applied to any point in the substrate and the resulting response determined. The cross-sectional mode is similar to the aerial mode except that the nodes are in a vertical cross section and gravity is included in the model. The substrate attachment is limited to the bottom and sides of the material region. The 3-D model is an extension into a third dimension of the 2-D model and deformation may be applied to the bottom and four sides of the material region.
The computational method used in solving for the deformation is a modified over-relaxation approach. In the basic over-relaxation method, equilibrium is reached by moving each node a distance proportional to the force acting on the node. The constant of proportionality is called the over-relaxation constant and optimal methods for selecting it are known. In each relaxation cycle, each node is moved and the relaxation cycles repeated until no node is moved further than a threshold distance, which is called the relaxation threshold. In the present invention, the over-relaxation is concentrated in those nodes where the greatest movement occurs. This significantly speeds up the computation time.
A novel aspect of the invention is the conditioning of the model. In the real world, the interpretation of the subsurface structure starts with observations of large scale fracturing and faulting. In one aspect of the invention, based upon the observed large-scale deformations, the user derives an initial geometry of the unfaulted material. Since it is clearly desirable that the end result of applying stresses to the model should at least duplicate the observed large-scale deformations, the model is weakened at the reconstructed fault positions so that upon application of the stresses, there will be a predisposition to produce the observed large-scale fractures and faults. The conditioning step is made easier in the present invention by use of the graphical user interface.
Another novel aspect of the invention is the use of an xe2x80x9canticipatexe2x80x9d step to get a quick solution to the deformation. The anticipate process is particularly useful in the definition of the deformation that is applied to the model. It is also useful in the conditioning of the initial trial model: comparison of the output of the anticipate step with the actual observed deformation serves as a check on the initial model. If the agreement is reasonable, then the initial trial model is used. If, however, the agreement is not good, a different initial trial model is used.